A new method for determining dipole-dipole energy in ID and 2D systems 
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An alternative method for computing dipole-dipole interaction energy in systems of ID and 2D 
periodicity like nanowires, nanotubes and thin films is presented. The approach is based on the use 
of periodic Green's functions that satisfy Laplace's equation and are analytically determined. The 
method, when combined with short-ranged interaction as in effective Hamiltonian, is suitable for 
studying finite-temperature properties of low-dimensional ferroelectric systems. 

PACS numbers: 



I. INTRODUCTION 

In the past decade it was shown that the finite- 
temperature behaviour of ferroelectric systems like 
BaTiC>3 can be successfully simulated by using the sta- 
tistical mechanics of an Effective Hamiltonian, which, in 
its turn, is based on the first-principles calculations of to- 
tal energies for small distortions of the high-temperature 
cubic structure. Such an approach predicts sequence of 
phase transformations, electromechanical responses and 
other finite-temperature properties in a good agreement 
with experiment. 

The conventional effective first-principles Hamiltonian 
approach is developed for the bulk systems infinitely re- 
peated in all three Cartesian directions (3D case) 0. 
Formally it can also be applied to the systems which are 
effectively infinite only in one or two dimensions sim- 
ply repeating them many times in finite direction(s) (in 
this case sufficiently thick vacuum gap(s) must be created 
between periodic replicas). Such a procedure, however, 
inevitably leads to errors in depolarizing electric fields 
and in corresponding shape-dependent electrostatic en- 
ergy [H 0. 0. 0| Consequently, in passing from the sys- 
tems with three- to that with two- and one-dimensional 
periodicity the part of the effective Hamiltonian con- 
nected with long-range dipole-dipole interactions should 
be modified. 

The existing methods for rapid evaluation of dipole- 
dipole interactions in ID 0,0 and 2D d 0, Eil EH s y s ~ 
terns are based on Ewald type summation technique ex- 
ploiting the integral representation of the Gamma func- 
tion (or Euler's integral) and Poisson summation for- 
mula. This technique leads to fast convergent sums in 
real and reciprocal spaces, which, however, are rather 
bulky. 

Here, we present a new elegant method for the treat- 
ment of dipole-dipole interaction in partially periodic sys- 
tems which leads to much simpler final expressions than 
the Ewald technique. The method is based on using the 
analytically determined periodic Green's function of the 
Laplace's equation Q(r',r). The simplicity of our final 
expressions is partly due to the fact that we do not use 
Ewald type transformations in ID case at all, whereas 
in 2D case we use it only for the special case when the 
the interacting dipoles lay in the same plane. The other 



reason comes from the fact that we perform the summa- 
tion of dipole-dipole interactions (from dipoles located in 
different unit cells) entirely in reciprocal space- with the 
only exception for the ID case when the dipoles lay along 
the line parallel to the periodicity direction. 



II. DIPOLE-DIPOLE INTERACTIONS IN 
SYSTEMS WITH ID AND 2D PERIODICITY 

We regard the systems with ID and 2D periodicity as 
infinite in one and two of the three Cartesian directions 
correspondingly. They can be obtained by replicating a 
" unit cell" an infinite number of times along that direc- 
tions. The dipole-dipole interaction energy can be writ- 
ten as 

1 x ^ f d(Ri) • d(Rj) 3 [d(R<) • Rij] [d(R,) • R, 



£dip — 9 X/ I 



Rf 



-J2J2 Da 0^ ~ R-jK(RiWRi), (!) 



a/3 i^j 



where = R; Rj, Rij = |Rjj | is the distance between 



dipoles, d a (Ri) is the a-component (a = x,y,z) of the 
dipole moment at the site i and 



D a p(Ri - Rj) = - lim 



_d d_( 



1 



o dr a drp \ \ 



R 7 ; + R, 



(2) 



Due to periodicity each vector R^ is represented as rj— R 



where is the dipole position inside the 0-th unit cell 
and R|| is an appropriate vector from the infinite number 
of vectors forming ID or 2D lattice. Accordingly, we may 
rewrite Eq.JI} as 



£dip = -£ Dapin-Tj+R^dairddpirj), (3) 

ck/3,R|| ij 

where the summation over i,j runs only inside the 0-th 
unit cell, the prime means that the term with = Vj in 
the case of Rm =0 must be omitted, N is the number of 
unit cells allowed to tend to infinity. We are interested 
in £ per unit cell, so N in Eq.© can be omitted. Using 
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this we can recast Eq.Q as 



1 ' d d 

£ dlp = - Jim - ^ ——G(r,r') | r ,= r( _ rj d a {r l )dp{r 3 ) 

a/3,ij a 

(4) 

where 



1 



R, ^ r - r ' + R W 



(5) 



is nothing but the periodic Green's function of Laplace's 
equation satisfying the point source equation 

V?0(r,r / )=-47r$^(r-r / + R||), (6) 

R i 



and the translation symmetry 

S(r,r') = S(r + R||,r' + R'||) 



(7) 



Now our task is to evaluate G(r, r') and then £di P for ID 
and 2D cases according to (@J. 



A. ID case 

Let z-axis to be along the infinite dimension, then all 
the vectors R|| lay in it. We represent each vector r as 
decomposed into two components {p, z}, where p is the 
projection of the r on the (x, y) plane. Accordingly, the 
Dirac functions <5(r — r' + R||) in the right-hand side of 
Eq.® become 5(p- p')S(z -z' + Ry). The 2D and the 
sum of ID Dirac functions can be expressed as Fourier 
integral and sum correspondingly: 



S(p-p') 



(27T) 



e *L-(p-p')dk. 



(8) 



£j(z-z' + R||) = I^ e lG ir( z - z '), (9) 



R,i 



where is the 2D wave-vector perpendicular to the z- 
direction, a is the period along this direction, and G|| 
are the reciprocal lattice vectors corresponding to the 
ID lattice of repeated cells. We shall look for the Green's 
function in the form 

S(r,r')=]T /. 9 (k ± ,G||) e t ( k -+ G ii>( r - r ')dk ± (10) 
G ii 

Inserting ©, ©, and (JTHJ into Eq.@ we obtainf!!!!^ 
g(r,r>) = - E / ~ 779 rr"? dk ± 



Gn + k^_ 



•]>>o(G|||p-p'|) 



(11) 



where Ko(x) is the 0-th order modified Bessel's function 
of the imaginary argument. Since the Bessel function 
Kq(x) decays exponentially for large argument x, the se- 
ries over G|| converge much faster than the direct 
lattice sum This is valid only for the contributions 
to Edip with p i pj. As to contributions coming from 
the dipole-dipole interactions with p i = pj or, in other 
words, from the chains of dipoles parallel to the z-axis, 
they can be easily calculated in the real space; the cor- 
responding ID sums are rapidly convergent like 1/n 3 . 
Moreover, doing so one can automatically solve the prob- 
lem of excluding of the "self-interaction" term from the 
sum ©. 

Substituting fTT|l into Eq.@J and taking into account 
that K' Q {x) = -K^x) and K'J{x) = K 2 {x) - Kx{x)/x 
one can easily obtain 

1 ' 

£di P = - e G \ cos ( G n ■ z y) 

x I K (G I pij ) d z (r i ) d z (i-j ) 

+ 77— Ki( G \\Pij) [d x (r l )d x (r j ) + ^(r^fo)] 

Lr|| Ptj 

-^ r K 2 {G\\p ij ) [p i3 • d(r,-)] [p i:j ■ d(r 3 -)] } 
Pij ) 

4EE g h sin ( G n • z «) K A G m) p-i? 

G,| ij 

[G,| • d(r,)] [ PtJ ■ d(r 3 -)] + [G|, • d(r 3 -)] [ Pij ■ d(r«)] 

^3 [dx (r< ) d x (r j ) + d y (Ti ) d v (r 3 ) - 2d z (r< ) d z (r 3 )] 

(12) 
Pj, Pi] = I Pij I, 



2a 3 
x/(Ny|) 



Here, G|| = |G|| 
and 



) °%3 



z i z j, Pij — Pi 



/(Kl)= E 



n— — 00 



"1.1 

a 



-3 



(13) 



Gn 



where the prime at the sum means that the term n = 
is to be excluded for the case i = j. The last sum in 
(|12fl describes the contribution to the energy associated 
with the chains parallel to the z-axis; this contribution is 
separated from the first two sums marked by the primes. 
Note that the function f(\zij\) is periodical in real space 
with period a. It is also worth noticing that the term 
C |l = does contribute to the sum i|12|) . Taking into 
consideration that as x — ► K (x) — + — ln(x), K±(x) — > 
1/x, and K%{x) — > 2/x 2 , we find 

£di P (G\\ = 0) = - E Pij 2 \ d x( r i) d x( r i) + d v( r i) d v( r j) 

ij 

-2p 4 r 2 [p y ..d(r j; )] [pydfo)]}. (14) 
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If all d(rj) || to z, this contribution turns to zero. In the 
macroscopic limit it describes the energy connected with 
the depolarizing field. 



B. 2D case 

Consider a slab or thin film with normal along the 
z-direction, each layer of which representing an infinite 
array of electric dipoles in the (x, y) plane. For this ge- 
ometry the vectors R|| form a 2D lattice parallel to this 
plane. The corresponding Green's function, as easily to 
show ac ting similar to the previous case and using the 
methods HE [yl, is 



e(r,r') = §]T 



2 ^ f e «k*-(z-z') e iG|| -(p-p') 



27T v-^ e -G||| Z - z '| e iG r (p-p') 

-E c n ' (15) 



Gil 



where S is the primitive unit sell area |15j . parallel to 
the slab plane, Gii are the reciprocal lattice vectors cor- 
responding to the 2D lattice and p = {x, y}. By inserting 
p5(l into Q it can be shown that only the contributions 
to Edip with Zi 7^ Zj exponentially converge as Gii in- 
creases (compare with the ID case). The terms with 
Zi = Zj describing the dipole-dipole interactions in the z- 
layers need additional care; we separate these terms and 
denote them as A. Now instead of Eq. (|I2|I we have 

E dip = 4 + -^ G \\ CX P(- G H Nl) 



Gii i,j 



B cos(G|| • Pij ) + C sin(G|| • Py)r J l 



. (16) 



where 



B = -d.Md.h) + [G,| • d(r,)] [G,| • d(r,)] , (17) 



C=- 



Gh 



[G r d(r i )]d,(r i )+[G r d(r J -)]d,(ri) . (18) 



The prime at the second sum in I jlfijl indicates the terms 
with Zij = to be dropped. To evaluate these terms 
(collected together in A) we apply the Ewald type trans- 
formation to the Green's function (|15|1 representing it in 
both coordinate and reciprocal spaces: 



ri >\ Y- erfc (^l r - r ' + R lll) . 20F 

Ry 

xVe iG l ■("-*') f V e-\ z - z '\ 2t2 e- G ii/ it2 t- 2 dt.(19) 



Here rj is the Ewald parameter, and erfc is the comple- 
mentary error function 



erfc(x) 



exp(-r) dt. 



(20) 



Using Ijl9|l and excluding the self-interaction term corre- 
sponding to Vi = Vj when Rm = 0, one obtains 

A = ?EE g h cos ( g h -pa) 



G,| ij 
[ 1 



1 



erfc 



G 



Git """ V ^ 
2T7 3 |d(r,)| 2 



G r d( ri )] [Gn-d^)] 



E 



30? 



(21) 



where the summation over i, j is constrained by the con- 
dition z^ — 0, r is the incomplete Gamma functional 



T(a, x) 



e-H a - x dt, 



(22) 



and the Ewald parameter -q is presumed to be large 
enough, so that the real space summation can be en- 
tirely neglected. It is interesting that the first sum in 
(|21ll contains the term Gm = 0. Indeed, since 



lim r --, 



1 Gl 



Gii ^0 



2' 4?y 2 



= 4r//G||, 



(23) 



we have 



A(G ll =0) = ^flY / d z (r i )d z (r j 

^Era 5 



(24) 



where D™ is the z-component of the total dipole moment 
in the n-th layer parallel to the (x, y) plane. 

Finally we note that in both ID and 2D cases the dipole 
energy can be rewritten as 

Edip = - E Qai3,i3d a (ri)dp(rj), (25) 

u(3,ij 

where Q a p,ij is the structure constant matrix; the latter 
can be calculated once and for all similar to the case of 
3D periodicity^]. Explicitly, the matrices are 
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1 1 2 

T K ^{ G Pii) P<*,ij PP,ij \ X] G Sill ( G ' K A G Pii) Pi/ X G aPfi,i] + 

P*3 > a G 

_^ oo t 

CI v 



n— — oo 

i 



1 G 2 



r I ~2' V J + G2 crfc ( ^ ) 



+G cxp(-G|zj 



,G a Gj} s . . G a 58 z . 

SazSpz) cos(G • p. — sin(G • p y ) — 



3v^F 



(26) 
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